{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create CMAQ-Ready Files from a Shapefiles\n",
    "\n",
    "This Notebook uses geopandas and cmaqsatproc to create IOAPI-like files for CMAQ. Geopandas supports optimized searches and projection conversions. This section focuses on creating I/O API NetCDF masks from **categorical** variables (e.g., **states**, **countries**) or from **quantitative** variables (i.e. **population**, **income**).\n",
    "* Categorical variables have grid cell values set to the fraction of the grid cell covered by a feature (e.g., state).\n",
    "* Quantitative variables have grid cell values set to the fraction of the feature (e.g., county) in the grid cell multiplied by the quantitative variable (e.g., population.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install Libraries\n",
    "\n",
    "- You can try uncommenting and then running this cell for these libraries to populate in your current kernel\n",
    "- If this doesn't work, you can try two different optons:\n",
    "  - try replacing !python -m pip` with `%pip`\n",
    "  - Or run pip install outside of the notebook for the python environment you will be using.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CWZSGMYnYHiD"
   },
   "outputs": [],
   "source": [
    "#!python -m pip install -qq cmaqsatproc geopandas xarray netcdf4 pycno"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libs Required\n",
    "- If this gives you a warning, your installation was likely unsuccessful\n",
    "- Or you need to restart your notebook to access the newly installed libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import shp2cmaq\n",
    "import warnings\n",
    "import cmaqsatproc as csp\n",
    "import geopandas as gpd\n",
    "import xarray\n",
    "import netCDF4\n",
    "import pycno\n",
    "# ignore warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0RQ1nVpN1-sL"
   },
   "source": [
    "## Create a CMAQ-Ready File from a Shapefile with Categorical Data (Part 1)\n",
    "\n",
    "This section uses shapefiles to create a CMAQ-ready file with variables for each feature (e.g., state or country). Each grid cell's value in a variable expresses the fraction of the grid cell area that is within the feature polygon (e.g., state or country)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download Shapefiles for Tutorial\n",
    "\n",
    "- You can try uncommenting the following !wget commands to pull in shapefiles\n",
    "- If the wget downloads don't work for you, try downloading the files from their respective websites:\n",
    "  - Natural Earth https://www.naturalearthdata.com.\n",
    "  - US Census https://census.gov/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DoAc3RR5B-vf"
   },
   "outputs": [],
   "source": [
    "# Download an example shapefile if you don't already have one.\n",
    "# default examples\n",
    "#!wget -N https://www2.census.gov/geo/tiger/GENZ2022/shp/cb_2022_us_state_500k.zip # use with attrkey STUSPS\n",
    "# alternate examples\n",
    "#!wget -N https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip # use with attrkey ADM0_A3\n",
    "#!wget -N https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip # us with attrkey iso_3166_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Configuration\n",
    "- Change the shppath to point to your shapefile\n",
    "- Set your attrkey (options for this should be in your shapefile metadata as a column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "M8SXK4c22Xwb",
    "outputId": "b46adad5-19c1-4330-ff65-8cf3a707a626"
   },
   "outputs": [],
   "source": [
    "# shppath : str\n",
    "#     Path to a shapefile or zip file containing a shapefile\n",
    "shppath = 'cb_2022_us_state_500k.zip'\n",
    "\n",
    "# attrkey : str\n",
    "#     Column to group shapes by e.g., STUSPS of census (AL, NC, etc)\n",
    "attrkey = 'STUSPS'\n",
    "\n",
    "# gdnam : str\n",
    "#     Name of grid definition within gdpath (e.g., 12US1, 108NHEMI2)\n",
    "gdnam = '12US1'\n",
    "gdpath = None # None uses built-in; or specify your own GRIDDESC path\n",
    "\n",
    "# For more options, run help(shp2cmaq.shp2cmaq)\n",
    "#?shp2cmaq.shp2cmaq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Tool\n",
    "\n",
    "1. This processor:\n",
    "    * Reads in native projection.\n",
    "    * Filters for in CMAQ domain.\n",
    "    * Optionally, custom extra processing. run help(shp2cmaq.shp2cmaq)\n",
    "2. Calculates area overlap:\n",
    "    * Performs grid cell intersections with shapefile polygons.\n",
    "    * Aggregates results to grid cell level.\n",
    "    * Finds largest area contributor.\n",
    "    * Calculates total cell overlap.\n",
    "3. Outputs:\n",
    "    * Store results as variables.\n",
    "    * Save as IOAPI-like file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "M8SXK4c22Xwb",
    "outputId": "b46adad5-19c1-4330-ff65-8cf3a707a626"
   },
   "outputs": [],
   "source": [
    "outpath = shp2cmaq.shp2cmaq(shppath, attrkey, gdnam, gdpath=None, verbose=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 383
    },
    "id": "i1wHo0m5s6j1",
    "outputId": "0173e1d0-0939-4d6e-ee96-f861b985bcb3"
   },
   "outputs": [],
   "source": [
    "igf = csp.open_ioapi(outpath)\n",
    "domkey = [k for k in list(igf.data_vars) if k.endswith('DOM')][0]\n",
    "totkey = [k for k in list(igf.data_vars) if k.endswith('TOT')][0]\n",
    "\n",
    "fig, axx = plt.subplots(1, 2, figsize=(12, 4))\n",
    "igf[domkey].where(lambda x: x > -999).plot(ax=axx[0], cmap='nipy_spectral')\n",
    "igf[totkey].plot(ax=axx[1], cmap='YlOrRd')\n",
    "_ = igf.csp.cno.drawcountries(ax=axx)\n",
    "name2idx = eval(igf[domkey].description)\n",
    "print('Dominant Index')\n",
    "print(str({v: k for k, v in name2idx.items()}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create CMAQ-Ready File from a Shapefile with Quantitative Data (Part 2)\n",
    "\n",
    "This section uses shapefiles to create a CMAQ-ready file with variables for each feature (e.g., state or county). Each grid cell in the variable expresses the fraction of a quantitative variable that is within that grid cell. The fraction within the cell is assumed proportional to the fraction of the feature polygon (e.g., state or county) in that grid cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download Shapefiles for Tutorial\n",
    "\n",
    "- You can try uncommenting the following !wget commands to pull in shapefiles and unzip them.\n",
    "- If the wget downloads don't work for you, try:\n",
    "    1. Go to https://censusreporter.org/,\n",
    "    2. Enter B01003 in the Explor dialog,\n",
    "    3. Enter states in the \"Show data by\",\n",
    "    4. Re-enter United States in the \"In\" dialog,\n",
    "    5. Hoover over \"Download data\",\n",
    "    6. Choose \"Shapefile\"\n",
    "    7. You must unzip the downloaded zip file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!wget -N -O acs2022_5yr_B01003_04000US21.zip 'https://api.censusreporter.org/1.0/data/download/acs2022_5yr?table_ids=B01003&geo_ids=040|01000US&format=shp'\n",
    "#!unzip acs2022_5yr_B01003_04000US21.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Configuration\n",
    "\n",
    "- Change the shppath to point to your shapefile **(for this you will first need to unzip your zip file and then point to your .shp file)**\n",
    "- Set your attrkey (options for this should be in your shapefile metadata as a column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shppath : str\n",
    "#     Path to a shapefile or zip file containing a shapefile\n",
    "shppath = 'acs2022_5yr_B01003_04000US21/acs2022_5yr_B01003_04000US21.shp'\n",
    "\n",
    "# attrkey : str\n",
    "#     Column to group shapes by e.g., STUSPS of census (AL, NC, etc)\n",
    "attrkey = 'geoid'\n",
    "\n",
    "# gdnam : str\n",
    "#     Name of grid definition within gdpath (e.g., 12US1, 108NHEMI2)\n",
    "gdnam = '36US3'\n",
    "gdpath = None # None uses built-in; or specify your own GRIDDESC path\n",
    "\n",
    "# srckey : str\n",
    "#     Name of quantitative variable (in this demo 'B01003001' is for population)\n",
    "srckey = 'B01003001' \n",
    "\n",
    "# For more options, run help(shp2cmaq.shp2cmaq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Tool\n",
    "\n",
    "1. This processor:\n",
    "    * Reads in native projection.\n",
    "    * Filters for in CMAQ domain.\n",
    "    * Optionally, custom extra processing. run help(shp2cmaq.shp2cmaq)\n",
    "2. Calculates area overlap:\n",
    "    * Performs grid cell intersections with shapefile polygons.\n",
    "    * Aggregates results to grid cell level.\n",
    "    * Finds largest area contributor.\n",
    "    * Calculates total cell overlap.\n",
    "3. Outputs:\n",
    "    * Store results as variables.\n",
    "    * Save as IOAPI-like file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "census_outpath = shp2cmaq.shp2cmaq(shppath, attrkey, gdnam, srckey=srckey)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "t3SJuEI0mFpV"
   },
   "outputs": [],
   "source": [
    "igf = csp.open_ioapi(census_outpath)\n",
    "domkey = [k for k in list(igf.data_vars) if k.endswith('DOM')][0]\n",
    "totkey = [k for k in list(igf.data_vars) if k.endswith('TOT')][0]\n",
    "\n",
    "fig, axx = plt.subplots(1, 2, figsize=(12, 4))\n",
    "igf[domkey].where(lambda x: x > -999).plot(ax=axx[0], cmap='nipy_spectral')\n",
    "igf[totkey].plot(ax=axx[1], cmap='YlOrRd')\n",
    "_ = igf.csp.cno.drawcountries(ax=axx)\n",
    "name2idx = eval(igf[domkey].description)\n",
    "print('Dominant Index')\n",
    "print(str({v: k for k, v in name2idx.items()}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3.6",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
